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We investigate the structure of a binary mixture of particles interacting via purely repulsive (point) 
Yukawa pair potentials with a common inverse screening length A. Using the hyper-netted chain 
closure to the Ornstein-Zernike equations, we find that for a system with 'ideal' (Berthelot mixing 
rule) pair potential parameters for the interaction between unlike species, the asymptotic decay of 
the total correlation functions crosses over from monotonic to damped oscillatory on increasing the 
fluid total density at fixed composition. This gives rise to a Kirkwood line in the phase diagram. 
We also consider a 'non-ideal' system, in which the Berthelot mixing rule is multiplied by a factor 
(1 + 5). For any 5 > the system exhibits fluid-fluid phase separation and remarkably the ultimate 
decay of the correlation functions is now monotonic for all (mixture) state points. Only in the limit 
of vanishing concentration of either species does one find oscillatory decay extending to r — oo. In 
the non-ideal case the simple random phase approximation provides a good description of the phase 
separation and the accompanying Lifshitz line. 

PACS numbers: 



I. INTRODUCTION 

A large class of fluids can be described generically 
as "big charged particles immersed in a neutralising 
medium of lighter particles" . Examples include charged 
colloidal suspensions 1 and dusty plasmas. 2 A simple 
model for such systems describes the effective interac- 
tion between the bigger particles in terms of a Yukawa 
(screened Coulomb) pair potential <fi(r) oc exp(— \r)/r. 
The effects of the screening due to the neutralising 
medium are incorporated via the screening parameter 
A. Such a Yukawa potential arises in, for example, 
the linearised Poisson-Boltzmann or Derjaguin-Landau- 
Verwey-Overbeek theories for the effective potential be- 
tween spherical charged colloids in solution. 1 The effect 
of the neutralising medium on the effective potential in- 
volves more than the screening effect, as described by 
the parameter A. There is an additional effect of charge 
renormalization, whereby the amplitude of the effective 
potential <fi(r) is not, as one might perhaps expect from 
a linear treatment, proportional to Z 2 , where Z is the 
charge on the big particles (colloids), rather the ampli- 
tude of <p(r) is proportional to Z 2 , where Z < Z is the 
renormalised charge i 

In the present paper we are concerned with a simple 
model of a binary mixture of big charged particles, with 
both species carrying charges of the same sign, immersed 
in a neutralising background medium. For a general bi- 
nary mixture of point Yukawa particles the pair poten- 
tials are dependent on six species specific parameters. 
These are the dimensionless coupling parameters, Mij, 
and the screening parameters, A^, for i,j = 1,2. We 
assume that both species experience the same screening, 
Ay = A, determined by the background medium (sol- 
vent), but that they have different coupling strengths. 
We can consider the coupling parameters to be pro- 
portional to the product of the effective charges of the 



species. Thus we define the pair potential as 



Myeexp(— Xr) 
AT 



(1) 



where e denotes the overall energy scale, and all three 
potentials are repulsive: My > 0. As usual it is as- 
sumed that the inter-species parameters are related to 
those for like particles. 3 The Berthelot mixing rule sets 
M\2 = \J M\\Mi2- This choice could correspond to ions 
with charge Zi, with Z%, Z2 > 0, immersed in a medium 
with an inverse screening length A. In order to generalise 
this mixing rule we introduce a non-ideality parameter 
6 > such that M 12 = (1 + (S)VMnM 22 . The case 6 = 
clearly corresponds to the ideal system. The non-ideal 
case <5^0 can be viewed as arising from charge screen- 
ing effects in the double layer of condensed counter ions 
on the surface of the particles. The double layer leads 
to an effective renormalisation of the particle charge and 
strongly affects the particle interactions^ One should ex- 
pect that for some positive (5 that the energy penalty in- 
curred for unlike species to be neighbours should lead to 
fluid- fluid demixing at high densities; such behaviour is 
found in models of soft-core fluids where positive non- 
additivity gives rise to demixing. 4 ' 5 ' 6,7 

Binary mixtures of H + and He ++ in a neutralising 
medium are expected to phase separate at temperatures 
and pressures of astrophysical interest - see Ref. 0. In 
this system the phase separation is thought to be due to 
charge neutralisation being less efficient in the mixture 
than in the pure phases. We believe that some of the 
complex screening effects associated with such systems 
may be incorporated in a simple model such as ours, via 
the parameter 6. Non-ideal charge renormalisation ef- 
fects may also be present in binary suspensions of col- 
loids. Charge renormalisation may be affected by the lo- 
cal concentrations of the different species of colloids and 
counterions in very subtle ways leading to the possibility 
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that the amplitude of <pi2 (r) may not be proportional to 
ZiZ 2 - Such effects would be mimicked using our simple 
model. It may be the case that in some colloidal fluid 
mixtures 6 > 0, or it may be that in other cases S < 0, 
i.e. non-ideal charge renormalization effects may favour 
mixing. However, in the present paper we investigate 
only the cases 8 = and S > 0. Using the accurate hyper- 
netted chain (HNC) approximation we find that for the 
case 6 = the total pair correlation functions hij (r) ex- 
hibit crossover from monotonic to exponentially damped 
oscillatory asymptotic decay, r — ► oo, on increasing the 
total density of the fluid mixture at fixed composition. 
This scenario is equivalent to that observed in the one- 
component point Yukawa fluid (OCY). 9 However, in the 
non-ideal case, an infinitesimal positive 6 can give rise 
to fluid-fluid phase separation at a sufficiently large to- 
tal density. Moreover we find that the fluid structure is 
changed profoundly from that pertaining to 5 = 0, i.e. for 
all thermodynamic state points, apart from the limits of 
pure species 1 and 2, the ultimate, r — > oo, decay of cor- 
relations is monotonic. We find that for the particular 
choice 8 = 0.1 the very simple random phase approxi- 
mation (RPA) provides a good account of the fluid-fluid 
binodal and spinodal obtained from the HNC but a poor 
account of the detailed behaviour of the correlation func- 
tions. 

The paper is arranged as follows: In Sees. HTland lllll we 
remind readers of the HNC and RPA integral equations 
and some basic results from the theory of the asymptotic 
decay of pair correlation functions in binary mixtures. 
Sec. IIVI describes results for hij (r) and for the poles of 
the Fourier transforms hij(q) obtained from numerical 
solutions of the HNC closure approximation. Within the 
RPA we are able to calculate the corresponding poles 
analytically and we compare these results with those from 
the HNC. The pole analysis enables us to determine the 
behaviour of hij(r) at intermediate range, as well as at 
longest range, r — * oo. Particular attention is paid to the 
behaviour of the correlation functions in the limit where 
the density of species 2, p 2 — > 0. In Sec. we present 
phase diagrams, in the composition, total density plane, 
along with Lifshitz lines for the partial structure factors, 
obtained from both the HNC and RPA. We draw some 
conclusions in Sec. IVII 



II. CLOSURE OF THE OZ EQUATIONS 

Our starting point for determining the fluid structure 
is the mixture Ornstein-Zernike (OZ) equation^ which 
relates the total correlation functions, hij (r) = gij (r) — 1, 
where gij{r) are the radial distribution functions, to a set 
of pair direct correlation functions, (r) : 



where Vij = |r$ — r 3 1 and pt is the bulk density of species 
k. These equations can be viewed as defining the pair 
direct correlation functions. In order to determine the 
fluid structure a second relation, or closure, is required. 

The simplest closure germane to the present model 
is the random phase approximation (RPA): c^j PA = 
—(3(j>ij(r) with j3 = {ksT)^ 1 , which is strictly valid 
only for r — > oo. Although this approximation is 
inadequate for hard-core model systems, it has been 
shown that the RPA becomes accurate for some soft- 
core systems at intermediate densities and exact at high 
densitiesA^SiLiSiiiiiS An important advantage of the 
RPA is that it does provide an analytical solution for 
correlation functions and for thermodynamic properties 
which may provide valuable physical insight into fluid 
behaviour. 

A more accurate approximation is the hyper-netted 
chain (HNC) approximation. The exact closure of the 
OZ equations can be expressed as: 

gij (r) = exp(-/30 i3 (r) + % (r) - Clj (r) - by (r)), (3) 

where — &y(r) is an unknown bridge function^ The HNC 
simply sets this bridge function to zero for all r. It 
is found to be accurate for long-ranged or soft-core 
potentials^2ii£ although it may fail in the neighbour- 
hood of a spinodal. For the one-component (point) 
Yukawa fluid it is known that the HNC is remarkably 
accurate for small coupling parameters In order to 
determine the correlation functions within the HNC we 
use a standard iterative procedure. In what follows, we 
fix the reduced temperature T* = (/3e) _1 so that the 
state of the system is determined by the total density, 
p = p± + p2, and the species concentrations Xi = Pi/p- 



III. ASYMPTOTIC DECAY OF CORRELATION 
FUNCTIONS 

There are two procedures for determining the asymp- 
totic, r — > oo, behaviour of the total correlation func- 
tions. One is to examine directly the numerical solutions 
for hij(r). The alternative method is to input the direct 
correlation functions (in the present case from either the 
RPA or HNC closures) into the set of of OZ equations @ 
and perform an asymptotic analysis. The OZ equations 
can be solved formally in Fourier space and the solution 
written as 



hij{q) 



D(q) : 



(4) 



htj(ri 2 ) = c i:j (n 2 ) + 



fc=l J 



d 3 r 3 c ifc (ri3)/i fc: ,(r32), (2) 



where hij (q) denotes the three-dimensional Fourier trans- 
form of hij(r). The three functions share the same de- 
nominator 

D (q) = [l - pica («)][! - P2c 2 2(q)} - piP2Ci 2 (q) 2 , (5) 
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but the numerators are dependent on the indices: 

Nu(q) = cn(q) + p 2 [ci2(q) 2 - cu(q)c22(q)], 
N 2 2{q) = 622(g) + Pi[ci2(q) 2 - C11 (5)622(5)], (6) 
N 12 (q) = N 21 (q) = c 12 (q). 

From the inverse Fourier transform it follows that: 



dqq sin(qr)hij(q). 



(7) 



Using Eq. (@J and assuming that the singularities of hij (q) 
for the present Yukawa systems are simple poles we are 
able to proceed via the residue theorem^ Performing 
contour integration around a semicircle in the upper half 
of the complex q plane, we write the total correlation 
functions as a sum of contributions from the poles en- 
closed: 



rhij (r) 



^2 A n exp(ig„r), 



(8) 



where q n satisfies D(q n ) = and A 1 ^ is the amplitude 
associated with the pole q n . This amplitude is related to 
the residue R% of qN ij {q)/D(q) by A% = R%/2tt. The 
poles are either purely imaginary, q = ia , or occur as a 
conjugate complex pair q = ±ai + iao^ 

In general there are an infinite number of poles and 
contributions from many of these are required to ac- 
count for the behaviour of hij(r) at small r. However, 
the ultimate, r — > 00, decay of hij{r) is determined 
by the pole that gives the slowest exponential decay, 
i.e. the pole with the smallest imaginary part. This is 
referred to as the leading order pole. If the leading or- 
der pole is purely imaginary then rhij(r) decays expo- 
nentially, rhij{r) ~ Aij exp(— Q-or), as r — > 00. On the 
other hand, if the leading order poles are a conjugate 
complex pair, then the sum of contributions from this 
pair of complex poles gives damped oscillatory decay, 
rhij(r) ~ 2Aij exp(— agr) cos(ctir — Oij), where Aij and 
Oij denote the amplitude and phase respectively^ 

Note that whereas the wavelength 2ir/ai and the de- 
cay lengths «q 1 or <Sq 1 are the same for all hij(r), the 
amplitudes and phases do depend on the indices ij*^ 
However, general considerations demand A\ 2 — ^11^22 
or A\ 2 = A11A22 and 28 12 = 6 n + «9 22 - 14 In the next 
section we shall employ the RPA and the HNC closures 
to investigate the ultimate, r — ► 00, decay of hij(r) and 
the behaviour of hij(r) at intermediate distances r. Note 
that the values of the amplitudes are relevant in deter- 
mining which pole or conjugate complex pair provides 
the dominant contribution to hij{r) in the intermediate 
regime. 



pA~ 3 = 0.5, x 2 — 0.5 and T* = 1. The coupling parame- 
ters are fixed at M\\ — 1, and M22 = 4. Fig. ^a) corre- 
sponds to ideal (Berthelot) mixing, (5 = 0, Fig. IHb) has 
5 = 10~ 5 , weak non-ideality, while Fig. ^c) has 6 = 0.1. 
In all cases the gij(r) appear structureless but the ex- 
panded scales of the insets reveal the nature of the inter- 
mediate and long range decay of the pair correlation func- 
tions. For <5 = 0, rhij(r) shows exponentially damped os- 
cillatory decay extending to arbitrarily large separations 
r. By contrast, for 5 = 10 -5 , there are several oscilla- 
tions at small r and rhij (f ) decays monotonically (expo- 
nentially) for Ar > 8. For <5 = 0.1 the monotonic decay 
extends from Ar ~ 3. Clearly the choice of 8 has a pro- 
found influence on the behaviour of the pair correlation 
functions. If one fixes the concentration at X2 = 0.5 and 
reduces the total density p one finds that for S = there 
is a crossover from the exponentially damped oscillatory 
behaviour of rhijir) shown in Fig.^a) to monotonic (ex- 
ponential) decay similar shown to that in Fig. ^c), near 
pA~ 3 ~ 0.05. Such crossover is found in the OCY on re- 
ducing the density at fixed T*»2> The scenario is different 
for (5 > 0. On reducing the density at fixed x 2 = 0.5 
the hij(r) are similar to those for pA~ 3 = 0.5; the decay 
remains monotonic. Similar features are found for other 
values of the concentration provided x±,x% ^ 0. If x% = 
or 1 then one recovers the OCY which exhibits crossover 
as mentioned above. 

In order to understand these results emanating from 
the full numerical solution of the HNC closure we turn to 
the asymptotic (pole) analysis of Sec. lTTTl It is convenient 
to begin with the simple RPA treatment before discussing 
the more sophisticated HNC results. 



A. Poles in the RPA 

The advantage of the RPA is that the pair direct corre- 
lation functions and their Fourier transforms, Cij(q), are 
given analytically. This means that the poles can be de- 
termined analytically. Using the definition of the RPA it 
follows from Eq. JD that 
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AT* 



(A 2 + (? 2 ) 



(9) 



Substituting this form into Eq. © we can solve for the 
zeros of D(q), i.e. the poles, q n , at each state point. We 
find that within the RPA there are only two poles. Both 
are purely imaginary and are given by 




M ± 



Ml 



M s 



A 2 . 



(10) 



IV. RESULTS FOR PAIR CORRELATION 
FUNCTIONS 

In Fig. ^ we display the radial distribution functions 
,j(r), obtained from the HNC closure, for the state point 



where we have introduced Mo = x\Mu + X2M22 and 
M s = Ax 1 x 2 M ll M 2 2{2 + 5)5. 

For <5 = this gives an imaginary pole with Q!q > 
A, pertaining to the positive sign, and a 'false' solution 
obtained with the negative sign giving ojq = A for all 
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FIG. 1: Pair correlation functions for the state point pA~ 3 
0.5. 12 = 0.5 and T* = 1, calculated from the HNC clo- 
sure, using the coupling parameter ratio M22/M11 = 4. Main 
figures show radial distribution functions, gij(r), and insets 
show In I r/iij(r) I versus Ar. a) S = 0. Exponentially damped 
oscillations extend to infinity, b) S = 1CF 5 . For Ar > 8, 
rhij(r) exhibit monotonic (exponential) decay, c) 5 = 0.1. 
Monotonic decay now develops for Ar > 3. 
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FIG. 2: The imaginary components of leading order poles, 
ao, calculated along a path in the phase diagram of increasing 
density, p, for fixed X2 = 0.5 and T* = 1. 5 = 0, correspond- 
ing to an 'ideal' mixture. The RPA solution, Eq. ill IH . con- 
sists of a single, purely imaginary pole, otg, which increases 
steadily with density from an /A = 1. Within the HNC clo- 
sure, at low densities p, the two leading order poles are both 
purely imaginary. As the density is increased the first imag- 
inary pole ascends and the second imaginary pole descends 
and, at the Kirkwood crossover point, coalesce becoming a 
conjugate complex pair q — ±«i + ian. Increasing the den- 
sity further both the real and imaginary parts of the complex 
poles increase. Thus the Kirkwood point marks the boundary 
between monotonic (at small pA -3 ) and damped oscillatory 
asymptotic decay. The inset shows the Kirkwood line, sepa- 
rating the two types of decay, plotted in the concentration- 
total density phase diagram. 



state points. The imaginary pole is analogous to that 
found in the RPA solution for the OCY. Fig.[21shows how 
this pole ascends from ao = A at p = as the density 
p is increased, along a line of constant concentration, 
X2 = 0.5, at fixed temperature T* = 1. 

By calculating the residues one can show (see Eq. (|14|) 
below) that the 'false' pole an = A makes no contribu- 
te, for 5 = the corresponding amplitude 



tion to hij (r) 



A,, = and 



rh« PA (r) 



A tj exp(-a+r) 



(11) 



where aj is given by Eq. with the positive sign, 

and the amplitudes are independent of the density^ 



A+= _MiL 
ij AT * 



(12) 



The situation is quite different for S > 0. We find that 
the pole with aj > A is modified slightly by the addition 
of a (relatively) small positive term. More significantly 
we find that the solution with the negative sign now cor- 
responds to a second pole which descends from olq = A 
eventually reaching zero as the density is increased - see 
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Fig. [3] The locus of points in the phase diagram for which 
Qq = is the spinodal, where the hij(q) (or the partial 
structure factors) diverge in the limit q — > 0, indicating 
that the system is undergoing fluid-fluid demixing. Note 
that this second 'spinodal' pole exists only for the mix- 
ture states, i.e. for Xi,X2 0; in the pure states it reverts 
to the 'false' solution ao = A. The behaviour of this pole 
is shown for fixed X2 — 0.5, and increasing density in 
Fig. El 

The 'spinodal' pole is present for arbitrarily small pos- 
itive 8. Therefore the RPA predicts that the system will 
undergo phase separation provided there is any degree 
of positive non- ideality in the pair potentials. This be- 
haviour is reflected in the asymptotic decay of the total 
correlation functions. For 6 > 0, 



rh 



RPA 

ij 



(r) = At exp(-o^r) + A tj exp(-a r) (13) 



with Q!q < A < ccq, which follows from Eq. (|10|l . The 
particularly simple forms of Eqs. I|llfl and are a con- 
sequence of the RPAii^ The amplitudes are 



2M 



iij (m ± y/M§ + M 4 ) + SijMg/a 
M Q + Ms/ (m ± y/M§ + M s 



A ± _ 1 
« 4AT* 



where superscript ± refers to the sign in front of the 
square root and we use the compound parameters Mq 
and Ms defined previously. Sij is the Kronecker delta. 
As the ultimate decay is determined by the pole with the 
smallest imaginary part, the 'spinodal' pole must de- 
termine the ultimate decay of hij(r) for any 5 > 0. For 
very small 8 > 0, lies very close to A until the reduced 
density pA -3 reaches very large values, i.e. the spinodal 
is shifted to very large densities as 8 — > 0. In these cir- 
cumstances the relative magnitude of the amplitudes A^- 
and A~j becomes important; note that a q depends only 
weakly on 8. For a given state point the amplitude of the 
'spinodal' pole A~- decreases relative to Afj as <5 — > and 
one must go to increasing values of r before the 'spin- 
odal' pole provides the dominant contribution to hijir). 
At intermediate r the contribution from aj determines 
the decay behaviour. 



B. Poles in the HNC 

For the OCY the HNC predicts a crossover line in the 
(p, T) plane separating a region where the asymptotic de- 
cay of h(r) is damped oscillatory from that where it is 
monotonic. The crossover takes place via the coalescence 
of two imaginary poles to form a complex pole as the 
density is increased at fixed temperature The mech- 
anism is equivalent to that discussed first by Kirkwood in 
pioneering studies of strong electrolytes^ Thus for the 
pure species we take the results of Ref. and read off the 
crossover values of pA~ 3 . For T* — 1 and X2 — (pure 
species 1) crossover occurs at pA~ 3 w 0.47 whereas for 



X2 = 1 (pure species 2) this occurs at pA -3 « 0.025. If 
we vary the concentration between these limiting values 
(at fixed T* = 1) we expect to find a crossover line in the 
(x2,pA -3 ) plane. We determined this line by calculating 
the poles of hij(q), i.e. the zeros of D(q), using the same 
numerical procedure as in Ref. © for the OCY. 

The pair direct correlation functions Cy(r) are ob- 
tained from the full solutions of the HNC integral equa- 
tions. In order to ensure convergence of the inte- 
grals which determine the poles we follow the proce- 



grals wmcn determine tlie po 
dure given in Refs. and ^3 



For r 



OO, Ci 



i(r) 



-Mij exp(— Xr)/(XT*r), so we define a short range di- 



rect correlation function c^(r): 



4 r (r) 



My exp(-Ar) 
? (r)+ XT*r ■ 



Fourier transforming we find 



eyte) = c-J(g) 



47tM,-. 



AT* (<z 2 + A 2 )' 



(15) 



(16) 



which can be substituted into Eq. JSJ. D{q) is separated 
into its real and imaginary components, and the equation 
D(q) = is solved numerically using a Newton- Raphson 
procedure. In general, the relevant integrals converge 
only for complex q such that ^s[q] < 2ao, where ao is the 
imaginary part of the leading order polei^i It follows that 
only a few poles can be calculated; the remaining poles 
are situated outside this region of convergence. Fortu- 
nately the poles relevant for determining the basic fea- 
tures of the decay of correlations can be obtained. 

For the case 8 = we find that for all X2 the leading 
poles exhibit Kirkwood cross-over mimicking that in the 
OCY. 9 At very low densities we find an imaginary pole 
(dashed line in Fig. |2J) that ascends from ao = A as we 
increase the density p at constant concentration, and is 
initially similar to the RPA. At slightly higher densities a 
second imaginary pole (dotted line) moves into the region 
of convergence and descends as p increases. These two 
poles move towards each other and then, at the Kirkwood 
point, coalesce. On increasing the density further, one 
finds a conjugate complex pair whose real and imaginary 
parts both increase with p (dash-dotted line). Fig. [21 
shows the imaginary parts of the relevant poles as they 
undergo the Kirkwood cross-over. The locus of points 
for which the poles coalesce, the Kirkwood line, is shown 
in the inset to Fig. [3 Below this line the asymptotic 
decay of r/iy(r) is pure exponential and above the line 
it is exponentially damped oscillatory. One sees that the 
state point X2 = 0.5 and pA -3 = 0.5, which corresponds 
to the results for hij(r) in Fig. QJa), lies deep within 
the oscillatory region so one expects to find oscillations 
extending to r — * oo. 

Fig. 01 displays the corresponding plots of ao for 8 = 
0.1. Now there are three purely imaginary poles at low 
densities. Two of these mimic closely what is found for 
(5 = 0, i.e, they coalesce and form a complex pole at a 
Kirkwood point that is not far removed from the corre- 
sponding point for (5 = 0. Once again one of these poles 
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FIG. 3: As in Fig. Hbut now for 8 = 0.1. Within the RPA 
the imaginary pole a g , is shifted by only a small amount from 
the case 8 = 0. A second purely imaginary pole is intro- 
duced which decreases from A to zero as p is increased; Oq = 
corresponds to the spinodal. Since < ctg this 'spinodal' 
pole determines the ultimate decay of the total correlation 
functions. Within the HNC closure the Kirkwood mechanism 
is still present and the crossover point is largely unchanged 
from the case 8 — 0. However, a 'spinodal' pole is intro- 
duced which follows closely the corresponding RPA pole . 
It follows that for all densities the ultimate decay of the total 
correlation functions is monotonic. Although the Kirkwood 
crossover is present it does not involve the leading order 'spin- 
odal' pole and therefore does not influence the ultimate decay. 



(dashed line) lies close to the RPA solution for very 
small values of pA~ 3 . The third pole follows closely the 
RPA 'spinodal' pole a^; it decreases slowly with increas- 
ing p\~ 3 for small densities before decreasing rapidly 
to zero at the spinodal point which is at pA~ 3 w 1.01, 
slightly higher than the value pA~ 3 = 0.98 found from 
the RPA. We refer to the third pole as the HNC 'spin- 
odal' pole. 

It is evident from Fig.Q|that within the HNC, as in the 
RPA, the purely imaginary 'spinodal' pole with ao < X 
will determine the ultimate decay of hij(r) which will be 
monotonic for all densities provided < x% < 1. The 
ultimate decay cannot be oscillatory. Although there is 
still a Kirkwood-likc crossover, this now involves higher 
order poles rather than the leading order poles, which was 
the case for 5 = 0. The corresponding crossover line is 
shown as the dotted line in the (x 2 , pA~ 3 ) phase diagram 
for S = 0.1 displayed in Fig.31 Note that the line lies far 
below the HNC spinodal (diamonds). Crossover does not 
influence the ultimate decay of hij (r) which is stubbornly 
monotonic. 

Of course, this does not mean that we do not find os- 
cillations at intermediate r. For densities larger than 
the crossover value the HNC pole analysis predicts, for 
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FIG. 4: Phase diagram for the binary Yukawa system with 
Mm/ Mn = 4, T* = 1 and 5 = 0.1. The system exhibits phase 
separation within both the RPA and HNC closures. The RPA 
spinodal (dashed curve) meets the RPA binodal (solid curve) 
at the critical point o. The HNC spinodal points (diamonds) 
were calculated by extrapolating the 'spinodal' pole to zero 
along lines of constant concentration and increasing density. 
The HNC binodal was calculated using the method outlined 
in the text; the straight segments are tie-lines connecting co- 
existing fluid phases denoted by crosses and correspond to the 
following reduced pressures (3P\~ 3 = 20, 22, 25, 30, 35, 40, 
45 and 50. Also shown is the Lifshitz line for Sjviv(g) cal- 
culated in both the RPA (dash-dotted) and HNC (crosses). 
The dotted line at low density denotes Kirkwood crossover 
from monotonic to oscillatory decay of correlation functions. 
However for 8 > and < xi < 1 the true asymptotic decay 
remains monotonic even for states above the line - see text. 

r — > oo, 

rhij(r) ~ Aij exp(— a^r) + 2A+j exp(— aor) cos(air — 0y) 

(17) 

with ao < okq. Provided Aij S> Aij and we are not close 
to the spinodal (ao = 0) the oscillatory contribution will 
be significant in the intermediate regime. 

An example is given in Fig. [S] which compares the full 
numerical result for /122 (r) with the contribution given by 
Eq. ijnjl for a statepoint (x 2 = 0.99, pA~ 3 = 0.05) where 
the amplitude of the 'spinodal' pole is very small. On 
a normal scale plot /122(f) is well-described by the con- 
tribution from a single conjugate complex pair of poles, 
i.e. by the second term in Eq. Q17[l. This contribution is 
oscillatory, albeit with a small amplitude. The inset to 
Fig.Oplots In |r/i22(f)| for a larger range of Xr. Now the 
oscillations are manifest but we observe monotonic (ex- 
ponential) decay, described accurately by the first term 
in Eq. (|T7Jl . for Xr > 14 (dotted line). As the ampli- 
tude XA22 of the 'spinodal' pole contribution is ~ 10~ 6 
this contribution does not become significant until large 
separations r. Of course such a contribution might be im- 
possible to detect experimentally or in simulations. We 
conclude that although a pair correlation function might 
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FIG. 5: Comparison between results of the full HNC solution 
for /122(f) (solid line) and the contributions of the lowest order 
poles for state point pA -3 = 0.05, X2 = 0.99 and parameters 
M22/M11 = 4, T* = 1 and 5 = 0.1. Note this point lies above 
the Kirkwood crossover line in Fig. 0] The main figure shows 
that the decay of /122(f) is well described by the contribution 
from a single pair of conjugate complex poles (dashed line) 
(ao/Ao = 2.269, Qi/A = 0.8723). On this scale the contribu- 
tion from the purely imaginary pole is vanishingly small. In 
the inset the plot of In |r/i22(r)| versus Ar (solid line) shows 
that although the oscillatory contribution (dashed line - vis- 
ible only in top left hand corner) describes /122(f) accurately 
for intermediate range, the contribution from the 'spinodaP 
pole (dotted line) (ao/A = 0.9993) determines the decay for 
Ar > 14. Thus the sum of contributions from the oscillatory 
and 'spinodal' poles (dot-dashed) given by Eq. I|170 is very 
accurate in the range 1.5 < Ar < 00. 



exhibit many oscillations, the ultimate decay can still 
be monotonic. Note that for a state point closer to the 
spinodal, such as in Fig.^c) where X2 = 0.5, p\~ 3 = 0.5 
and ao/A ~ 0.72, the purely exponential contribution in 
Eq. Ijl7|l dominates for Ar > 3. 

For a pure (one component) fluid there is no 'spinodal' 
pole; within the HNC the Kirkwood crossover occurs be- 
tween leading-order poles. Therefore, on the axes X2 = 
and xi = in Fig. ^hu(r) must decay in an oscillatory 
fashion as r — > 00, provided we consider states above the 
Kirkwood crossover densities. However, as we increase X2 
(or x\) infinitesimally the 'spinodal' pole makes a non- 
vanishing contribution to the correlation functions and 
forces the ultimate decay of hij{r) to be monotonic. In 
order to understand the evolution of this behaviour we 
must consider the amplitudes, A^, of the 'spinodal' pole 
as the concentration X2 vanishes along a line of constant 
total density. Fig. EJshows results for p\~ 3 = 0.5, a den- 
sity that lies just above the Kirkwood crossover value, 
p\~ 3 ~ 0.47, for the pure fluid. The 'spinodal' pole 
c*o/A — > 1~ as X2 — ► 0, and the accompanying amplitudes 
An (positive) and A12 (negative) both vanish, as power 
laws, in this limit. It follows that the 'spinodal' pole 
makes contributions to hu(r) and /112(f) that become 



FIG. 6: The HNC 'spinodal' pole qo and corresponding am- 
plitudes Aij as a function of concentration X2 along a path 
of constant density pA~ 3 = 0.5, for T* = 1, and 5 — 0.1. As 
x 2 -> 0, a /A — > 1", Aii,^^ — > but A22A tends to 0.84. 
Note that A12 is negative and its magnitude is plotted here. 
The amplitudes obey the rule Ai 2 = A11A22 - see text. 



vanishingly small in a continuous fashion, as a; 2 — > 0. 
Thus for very small concentrations hn(r) and /112(f) ap- 
pear oscillatory until extremely large values of Ar where 
the 'spinodal' pole will dictate the final (monotonic) de- 
cay of these functions. 

The amplitude A22 of the minority component, species 
2, has a different variation. As shown in Fig. EI A^22 re- 
mains almost constant as the concentration is reduced 
and tends to a non-zero value (0.84) at X2 — 0. At first 
sight this result is a little surprising. However we should 
recall that physical observables such as probability distri- 
butions or the liquid structure factors involve the product 
P2^22(f) which vanishes as X2 — ► 0. We confirmed that 
our numerical results for the amplitudes satisfy the rule 1 ^ 
A\ 2 = A11A22 mentioned in Sec. IIIII 

At this point it is instructive to return to the RPA 
'spinodal' pole. Recall that is given by Eq. (|10(l and 
the corresponding amplitudes by Eq. (|14f> . If we Taylor 
expand around X2 = we find that 



= A 1- 



2n P M 22 5(2 + 6) 
X 3 T* 



x 2 )+0(x A 2 ), (18) 



j4 1:l and A 12 decay to zero as power laws in X2, and that 



A 22 tends to a constant, i.e. 



A11 = 



Ml 2 (l + S) 2 (2 + S)S 
AT*Mn 



xl + 0{x%), 



Ml 2 /2 {\ + 8){2 + 5)5 



1 12 



i 22 



AT*M 1 1 1 /2 



x 2 +0(x 2 2 ), (19) 



In the limit X2 



M 22 (2 + (?)<? 
AT* 

0, XAZ 



0(x 2 ) 



M 22 (2 + 5)8 /T* which 



takes the value 0.84 for the present choice of parameters: 
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M 2 2 = 4, S = 0.1 and T* = 1. This is the same limiting 
value of the amplitudes as obtained in the HNC. Indeed 
we find that the RPA results for and for the amplitude 

A~- are almost identical to those from the numerical HNC 

1 j 

results for concentrations x 2 % 10~ 4 . Thus, numerically, 
the limiting behaviour predicted by Eqs. (fTH|l and ijT§|) 
pertains to the HNC as well as to the RPA. Note that 
the coefficients of the leading order terms in Eq. I(19|) are 
consistent with the requirement (A^ 2 ) 2 — A^ X A^ 2 . 

It is not immediately obvious why the two (very) differ- 
ent closure approximations should yield the same limiting 
behaviour for the pair correlation functions as x 2 — > 0. 
That they do suggests that the limiting behaviour de- 
scribed by Eqs. (|18|) and l|19l) should be valid more gen- 
erally. In the Appendix we show that these results fol- 
low from the natural division, Eq. 1)15(1 . of the pair di- 
rect correlation functions into a Yukawa tail plus a short 
ranged contribution cfj(r). Provided the Fourier trans- 
forms cfj(g) are sufficiently well-behaved in the neigh- 
bourhood of q = iX, Eqs. (|18fl and (|19|) should be valid 
for any closure. 



V. PHASE DIAGRAM AND LIFSHITZ LINE 

FOR 8 = 0.1 

The results presented in the last section imply that 
within the RPA and HNC introducing a small degree of 
non-ideality into the pair potentials defining the mixture, 
i.e. imposing 6 > 0, leads to a 'spinodaP pole which dic- 
tates that pair correlation functions should decay mono- 
tonically as r — > oo for all state points, including those 
far removed from the spinodal, where one might have 
expected oscillatory decay (within the HNC). In this sec- 
tion we focus on the spinodal itself and the associated 
fluid-fluid phase separation. 

The spinodal is the locus of points in the phase dia- 
gram at which the 'spinodal' pole reaches zero. Within 
the RPA we can determine the spinodal by finding so- 
lutions of <2q =0, where a$ is given by Eq. (fTU|> with 
the negative sign. This result is plotted in Fig. El for 
8 = 0.1 and T* = 1 (dashed curve). Within the HNC it 
is not possible to determine the spinodal precisely. For 
certain densities and concentrations the HNC approxi- 
mation does not have solutions. However it is possible 
to calculate the position of the HNC 'spinodal' pole as 
a function of p\~ 3 at fixed concentration and extrapo- 
late to zero as shown in Fig. [3] These spinodal points, 
obtained for a number of concentrations, are shown in 
Fig- El (diamonds). The HNC spinodal lies close to that 
obtained in the RPA. 

Using the RPA we are able to write the reduced bulk 
Helmholtz free energy per particle, /, as a sum of ideal 
and excess parts^Sii 

f(p,x 2 ) = fid + ^p[xUn(0) + 2a; 1 a; 2 ^ 12 (0) + xj$ 22 {0)] 

(20) 



where <f>ij(0) — ineMijX^ 3 is the q = limit of the 
Fourier transformed pair potential. The ideal part, 
fid(p,%2), contains the ideal free energy of mixing, 
f3~ 1 {xi ln(xi) + x 2 ln(a; 2 )}, as well as a term in p that 
is irrelevant for phase behaviour. Eq. (|20l) corresponds 
to calculating / from the compressibility route, and the 
spinodal obtained from this equation is identical to that 
obtained from the zeros of . We now Legendre trans- 
form to the Gibbs free energy g = f + Pv, where v = 1/p 
is the volume per particle and the pressure is given by 
P = —(df/dv) X2 . Using the common tangent construc- 
tion on g we obtain the binodal which is also plotted in 
Fig- El (solid curve). The binodal and spinodal meet at 
the critical point near x 2 = 0.27, pA~ 3 = 0.88. 

In implementing the HNC we choose to calculate the 
thermodynamic functions locally, thereby avoiding ther- 
modynamic integration. We begin by tracing out isobars 
across the phase diagram with the pressure given by the 
virial route: 

i j 

(21) 

Along the isobars we calculate the HNC chemical poten- 
tials 

PlH = ln(Mi) (22) 
+ ]T Pj J dr | [hij (r) - dj (r)] - c« (r) } 

where Aj is the (irrelevant) thermal de Broglie wave- 
length of species i, with i = 1,2. The Gibbs free energy 
per particle is then given by g — x\p\ + £ 2 /i 2 . By fit- 
ting a polynomial to these results we are then able to 
determine the binodal using the common tangent con- 
struction. For the isobars that intersect the region for 
which the HNC does not provide a solution we calculate 
the two 'branches' of the free energy and construct the 
common tangent on these; the procedure is equivalent to 
that used in Ref. 0- Nevertheless, for state points close 
to the critical point it is not possible to determine the 
binodal. The HNC binodal is shown in Fig. El (crosses 
connected by tie lines). This lies outside our estimated 
HNC spinodal and fairly close to the RPA binodal. A 
slightly better level of agreement between the HNC and 
RPA results for the binodal was found in Ref. |a for a 
softcore model of binary star-polymers. 

Since the 'spinodal' pole dictates the ultimate decay of 
the pair correlations even for states that are far removed 
from the spinodal it is instructive to seek some crite- 
rion which indicates when contributions from the other 
pole(s) become important in determining the structure 
of the fluid. One valuable indicator of the change in the 
latter is the Lifshitz line which focuses on the behaviour 
of the fluid structure at small wave numbers qfi ,l $ We 
define the partial structure factors asi^ 

Sij(q) = S. Lj + (xiXj) 1/2 phij(q), (23) 
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and concentrate on the number-number structure factor 

S NN (q) = x 1 S 11 (q)+x 2 S 22 (q) + 2(x 1 X 2 ) 1 / 2 S 12 (q). (24) 

As we approach the critical point or, indeed, the spin- 
odal, the partial structure factors diverge at small q: 
Su(q = 0) and S 22 (q = 0) +oo and S 12 (q = 0) 
— oo. These divergences are reflected in the linear com- 
bination, i.e. SpfN^q = [))—► +co. 

We define the Lifshitz line as separating regions of the 
phase diagram for which S^Niq) has a local maximum 
or a local minimum at q — 0. We choose Snn(i) because 
i) it is the most symmetrical combination of the partial 
structure factors, and ii) it diverges in a similar manner 
on both sides of the phase diagram. Formally we make a 
small q expansion of Snn{q)' 

S NN {q) = a(p, x 2 ) + b(p, x 2 )q 2 + 0{q A ), (25) 

so that the Lifshitz line is the locus of points for which 
b(p,x) = 0, 6 ^ This line is calculated in both the RPA 
and HNC and the results shown as the dash-dotted curve 
and crosses, respectively, in Fig. We find that there 
is good agreement between the two closures, suggest- 
ing that in this region of the phase diagram the RPA 
should become a fairly reliable approximation for describ- 
ing hij(r) at both intermediate and large values of r. 

VI. DISCUSSION AND CONCLUSIONS 

In this paper we have investigated the structural and 
thermodynamic properties of binary mixtures of parti- 
cles interacting via purely repulsive (point) Yukawa pair 
potentials. We have used two approximate closures to 
the OZ equations in order to calculate the fluid proper- 
ties: the simple RPA allows us to elucidate analytically 
some important properties, and the HNC approximation 
is expected to be very accurate for the Yukawa fluid at 
the intermediate fluid densities relevant to the present 
work. This expectation stems from the fact that the HNC 
closure correlation functions are rather accurate for the 

For the 'ideal' mixture, with 5 = 0, the RPA approx- 
imation gives a poor description of the fluid structure, 
predicting monotonic (exponential) decay of r/iy (r) for 
all state points in the phase diagram rather than the 
Kirkwood crossover line that is manifest in the HNC - 
see Fig. |2 No fluid-fluid phase separation occurs for 
6 = 0. By contrast, when 8 > an additional purely 
imaginary pole arises within the HNC approximation for 
hij(q). This 'spinodal' pole governs the phase separa- 
tion in the mixture for the non-ideal case. Since the 
RPA provides a good account of this pole - see Fig. 0- 
it follows that the RPA also yields fluid-fluid transition 
boundaries that are quite close to those from the more 
accurate HNC - see Fig. The amplitude Aj- of the 
spinodal pole contribution to h%j (r) becomes increasingly 
large with increasing proximity to the spinodal, and this 



pole provides the dominant contribution in the region of 
the phase diagram near the spinodal (roughly the region 
enclosed by the Lifshitz line - see Fig.^J. 

For states away from the spinodal/binodal, particu- 
larly in the limits x\ — > or x 2 — > 0, the amplitudes of 
the contribution from the spinodal pole to the majority 
species correlation function hu(r) (where i is the major- 
ity species) and to hi 2 (r), become very small - see Figs.O 
and El One must magnify the large r tail of these corre- 
lation functions in order to see the contribution from the 
spinodal pole. Above the Kirkwood line, but well away 
from the spinodal, the intermediate r decay of hij(r) is 
damped oscillatory; the oscillations can persist to large 
r, before the monotonic decay from the spinodal pole fi- 
nally dominates at very large values of r. This result 
is quite unusual, since one generally expects the leading 
order pole to dominate the decay of (r) at both large 
and intermediate values of r^LIiiil In the limit x 2 — ► 0, 
the spinodal pole — > iX and the amplitudes, and 
Ai 2 , of the contribution of to h\\{r) and to hi 2 (r) have 
a power law dependence on x 2 (see Eq. (|19|l L vanishing 
in the pure fluid of species 1. However, the amplitude, 
A^ 2 , of the contribution from the spinodal pole to h 22 (r) 
tends to a non-vanishing value in the limit x 2 — > 0. This 
means that for small values of x 2 the nature of the decay 
of the correlation function h 22 {r) can be very different 
from that of hn(r) and h\ 2 {r), a situation which has not 
been encountered in other fluid mixtures. In experimen- 
tal or simulation results it would be almost impossible 
to detect the contribution in hi\{r) and hi 2 (r) from the 
'spinodal' pole for small x 2 . However, its existence could 
be inferred from the apparently different decay of h 22 (r) 
versus that of hn(r) and h\ 2 (r). 

It is important to emphasise that the amplitude A^ 2 
tending to a non-zero value in the limit x 2 — > does 
not imply any pathological consequences for thermody- 
namic or measurable structural quantities. As pointed 
out earlier, physical observables such as the liquid struc- 
ture factors involve the product p 2 h 22 (r) which vanishes 
as a; 2 — > 0. One exception to this rule is the effective 
potential between two particles immersed in a solvent of 
the other species. The effective potential 022^ ( r ) between 
two particles of species 2 at infinite dilution is related to 
the radial distribution function, in the limit x 2 — » 0, by 

g 22 (r)=eM-^t f {r))- (26) 

Clearly (jx^J (r) depends on the density and temperature 
of the solvent (species 1). We define the solvent- mediated 
potential W 22 (r) via $ 22 if) = 4> 22 (r) + W 22 (r), where 
4> 22 {r) is the bare (direct) potential between the two par- 
ticles. W 22 (r) depends on the nature of the solvent and 
on the solvent-particle interaction. From Eq. (12(:i[l it fol- 
lows that W 22 (r) depends on h 22 (r) rather than p 2 h 22 (r). 
There are clear implications for the form of 022^ ( r ) given 
that the amplitude of the contribution to h 22 (r) from the 
spinodal pole is non-vanishing in the limit x 2 — > 0. We 
leave a full discussion to another publication but point 



10 



out that the presence of the spinodal pole, which occurs 
for any S > 0, gives rise to an effective potential that 
is attractive at sufficiently large distances and decays as 
cxp(— Xr)/r, i.e. with the same length scale as the bare 
potentials. 

Given that the Yukawa potential arises in a great vari- 
ety of physical problems, si typically where there are big 
charged particles screened by a neutralising background 
medium, we believe the present work should be relevant 
to systems ranging from charged colloidal fluids! to dusty 
plasmas, 2 

We have chosen to focus on a simple example of a bi- 
nary Yukawa mixture, namely one in which the inverse 
screening length A is the same for all ij pairs, as this is the 
situation which arises naturally for charged particles im- 
mersed in a neutralising medium. One could, of course, 
consider more complex examples with A dependent on 
the species indices Moreover, one could consider mix- 
tures that correspond to negative non-additivity, i.e. with 
5 < 0. In this case Mn < \/M\\M22 and the interactions 
should favour fluid fluid mixing. 
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Appendix: Behaviour of correlation functions as 

x 2 -> 

Here we show that the results for the spinodal pole and 
the amplitudes AZ, Eqs. l(T%)l and ifHfl) . in the limit x-i — > 
are not specific to the RPA but follow generally from 
the form of the direct correlation functions for Yukawa 
mixtures. Making the separation of the Cy (r) given by 
Eq. i|15|) and substituting into Eq. JSJ we obtain 



and ctij = AirMij / XT* . The poles are given by the solu- 
tion to the equation D(q) = 0. From equation 1)27(1 . we 
see that one set of solutions is given by 



p± = -(b±\/b 2 - 4ac)/2a. 
a and b are functions of q. However, Eq 
purely imaginary poles at q± = iao = iyX? 



(29) 

(121 leads to 



p±, pro- 



vided we assume that the functions cfj(g) are well be- 
haved (finite and differentiable) on the imaginary axis 
around q± . The leading order pole is that corresponding 
to p^, which in the limit of vanishing density p2 = 
(i.e. c = 0) yields a pole at q = iX, giving a decay 
rhij(r) ~ A~j exp(— Ar), where the amplitudes A~- are to 
be determined below. For small concentrations of species 
2 we can Taylor expand p_ in powers of c, giving 



P- = 



piotxx 



c(b - giOu) 

2 2 

Pi^ii 



ac 
Pi"li 



O(pl). (30) 



We find that in the limit p2 — > 0, p- ~ —(0:22 — 
a i2/ a ii)/°2- Thus the leading order pole has 



which is identical to the RPA result in Eq. (|18fl . 
The amplitude, A^- is given by: 



a: 



q-Nijjg-) 
2vrL>'(g_) ' 



(32) 



where the prime denotes the derivative with respect to 
q and the functions Nij(q) are given in Eq. ©. From 
Eq. J23) we find that 



D(q) 



b c 

aH h - 

P P 



(27) 



D'(q) = a' + 



P 



2qb 

P 2 



Aqc 



P 



(33) 



where 

P 
a 
b 



q 2 + X\ 



[1 - pxc\\{q)][l - p 2 c s 2 r 2 (q)] ~ PxP2K 2 (q)f 
[1 - pxc s 1 r 1 (q)}p 2 a 2 2 + [1 - p2C 22 {q)]pia 11 
+2p 1 p 2 c s 1 r 2 (q)a 12 , 
PiP2(aiia 2 2 - a\ 2 ) 



In the limit P2 — > 0, using Q30[l . we find that D'(q_) = 
2q r _&/p 2 _+other terms less singular in p2- Using this re- 
sult, we can expand in Eq. 132|) to obtain the amplitudes 
A~ in the limit P2 — > 0. We find that the leading order 
terms are precisely the RPA results given in Eq. I|19|) - 
i.e. in the limit P2 — > the amplitudes A~- are indepen- 
(28) dent of the functions cfUq). 
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